Björn Reineking, 2022-09-09
Model the relative density of animals (also called range distribution or utilisation distribution) as a function of environmental predictors.
We will use the buffalo data set. Loading packages
library(animove)
library(ctmm)
library(sf)
library(mgcv)
library(mvtnorm)
library(lubridate)
See Kami’s slides for how we got here…
data(buffalo_utm)
data(buffalo_env)
raster::plot(buffalo_env)
raster::plot(raster(buffalo_env, 1))
points(buffalo_utm)
To speed up analyses, we will only work with one individual, Cilla
cilla <- buffalo_utm[["Cilla"]]
raster::plot(raster(buffalo_env, 1), ext = extent(cilla) * 2)
lines(cilla)
The first day of Cilla shows some unrealistic movements. We remove those observations
cilla <- cilla[timestamps(cilla) > min(timestamps(cilla)) + days(1), ]
raster::plot(raster(buffalo_env, 1), ext = extent(cilla) * 2)
lines(cilla)
Create telemetry object
cilla_telemetry <- as.telemetry(cilla)
## Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
## to PROJ string
## Warning in wkt(obj): CRS object has no comment
## Warning in wkt(obj): CRS object has no comment
## Minimum sampling interval of 3 minutes in Cilla
Fit ctmm model
cilla_guess <- ctmm.guess(cilla_telemetry, CTMM=ctmm(isotropic = TRUE), interactive = FALSE)
cilla_select <- ctmm.select(cilla_telemetry, cilla_guess)
Fit akde
cilla_akde <- akde(cilla_telemetry, cilla_select)
plot(cilla_akde)
Create named list of rasters We can do this by hand
be <- list("elev" = raster(buffalo_env, "elev"),
"slope" = raster(buffalo_env, "slope"),
"var_NDVI" = raster(buffalo_env, "var_NDVI"))
The integrator = “Riemann” option is still experimental, we use it here because it is much faster
cilla_rsf_riemann <- rsf.fit(cilla_telemetry, cilla_akde, R = be, integrator = "Riemann")
## Warning in .local(x, ...): This function is only useful for Raster* objects with
## a longitude/latitude coordinates
## Maximizing likelihood @ n=166037
## Calculating Hessian
summary(cilla_rsf_riemann)
## $name
## [1] "OUF"
##
## $DOF
## mean area diffusion speed
## 8.322854 16.577732 17.688692 17.345680
##
## $CI
## low est high
## var_NDVI (1/var_NDVI) -188.04892994 -72.68984716 42.66923562
## slope (1/slope) -0.37986892 -0.03499126 0.30988641
## elev (1/elev) -0.04693589 -0.01277037 0.02139514
## area (square kilometers) 230.24198135 398.32893545 611.73047276
## τ[position] (days) 4.56389303 7.82598060 13.41967744
## τ[velocity] (minutes) 40.67981533 43.18035421 45.83459818
## speed (kilometers/day) 10.35589660 13.52672089 16.69123828
## diffusion (square kilometers/day) 3.12109675 5.29364725 8.03075027
A suitability map
suitability_riemann <- suitability(cilla_rsf_riemann, be, crop(be[[1]], extent(cilla) * 2))
raster::plot(suitability_riemann)
Range distribution (includes the ranging behaviour)
agde_cilla <- agde(cilla_rsf_riemann, be)
plot(agde_cilla)
rsf.predict <- function(model, object, include_avail = TRUE) {
log_avail <- 0
if(include_avail) {
xy <- as.data.frame(xyFromCell(object, 1:ncell(object)))
xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(object))
xy <- st_transform(xy, CRS(model@info$projection))
xy <- st_coordinates(xy)
xy <- cbind(xy, -(xy[,1]^2 + xy[,2]^2)/2)
beta_rr <- 1/model$sigma[1,1]
beta_xyr <- c(model$mu*beta_rr, beta_rr)
log_avail <- xy %*% beta_xyr
}
newdata <- as.data.frame(getValues(object))
attr(newdata, "na.action") <- "na.pass"
X <- stats::model.matrix.default(model$formula, data = newdata, na.action = na.pass)
lambda <- X[,-1] %*% model$beta + log_avail # model.matrix.default puts as first column (Intercept) if model formula does not include "-1"
r <- raster(object, 1)
r[] <- exp(lambda - max(lambda, na.rm = TRUE))
r <- r / sum(getValues(r), na.rm = TRUE)
r
}
cilla_predict <- rsf.predict(cilla_rsf_riemann, crop(buffalo_env, extent(cilla) * 2), include_avail = FALSE)
raster::plot(cilla_predict)
Functions to generate quadrature points and predict with the model
rsf_points <- function(x, UD, R = NULL, n = 1e5, k = 1e6, type = "Riemann",
rmax =6*sqrt(UD@CTMM$sigma[1,1]),
interpolation = FALSE) {
# Samples background points from a 2D normal distribution fitted to the relocation data, and extracts environmental
# information from a raster object "R"
# x: telemetry object
# UD: UD object
# R: raster* object
# n: number of background points to sample
# k: weight of presence points
# rmax: maximum distance for Riemann-type integration
# interpolation: do interpolation when sampling the grid
# When type 0´= "MonteCarlo", importance sampling is done
stopifnot(UD@CTMM$isotropic)
stopifnot(type %in% c("Riemann", "MonteCarlo"))
if (type == "Riemann") {
quadrature_pts <- getValues(R)
xy <- as.data.frame(xyFromCell(R, 1:ncell(R)))
xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(R))
xy <- st_transform(xy, CRS(UD@info$projection))
xy <- st_coordinates(xy)
r <- sqrt(((xy[,1] - UD@CTMM$mu[1]))^2 + ((xy[,2] - UD@CTMM$mu[2]))^2)
bg <- data.frame(case_ = 0,
x_ = xy[r<rmax,1], y_ = xy[r<rmax,2], w_ = prod(res(R)), k_ = k)
bg <- cbind(bg, quadrature_pts[r<rmax,])
bg <- sf::st_as_sf(bg, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
xx <- data.frame(case_ = 1, x_ = x$x, y_ = x$y,
w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k)
xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
xx[names(R)] <- as.data.frame(raster::extract(R, st_transform(xx, crs(R)), method = ifelse(interpolation, "bilinear", "simple")))
xx <- rbind(bg, xx)
} else {
quadrature_pts <- MASS::mvrnorm(n, mu = UD@CTMM$mu, Sigma = UD@CTMM$sigma)
xx <- data.frame(case_ = 0, x_ = quadrature_pts[, 1], y_ = quadrature_pts[, 2], w_ = UD@CTMM$sigma[1,1]/n, k_ = k)
xx <- rbind(xx, data.frame(case_ = 1, x_ = x$x, y_ = x$y,
w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k
))
xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = CRS(UD@info$projection))
xx[names(R)] <- as.data.frame(raster::extract(R, st_transform(xx, crs(R)), method = ifelse(interpolation, "bilinear", "simple")))
}
xy <- st_coordinates(xx)
colnames(xy) <- c("x_", "y_")
sd <- sqrt(UD@CTMM$sigma[1,1])
xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
xx <- cbind(xy, xx)
xx
}
predict_rsf <- function(model, UD, object, include_avail = TRUE, data_crs = UD@info$projection) {
if(include_avail) {
xy <- as.data.frame(xyFromCell(object, 1:ncell(object)))
xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(object))
xy <- st_transform(xy, data_crs)
xy <- st_coordinates(xy)
colnames(xy) <- c("x_", "y_")
sd <- sqrt(UD@CTMM$sigma[1,1])
xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
} else {
xy <- cbind("x_" = rep(0, ncell(object)), "y_" = rep(0, ncell(object)))
}
newdata <- as.data.frame(cbind(xy, getValues(object)))
lambda <- as.numeric(predict(model, newdata, type = "link"))
r <- raster(object, 1)
r[] <- exp(lambda - max(lambda, na.rm = TRUE))
r <- r / sum(getValues(r), na.rm = TRUE)
r
}
Generate quadrature points (“background points”)
set.seed(2)
rsf_cilla_df <- rsf_points(cilla_telemetry, cilla_akde, buffalo_env, interpolation = TRUE)
rsf_cilla_df <- rsf_cilla_df[!is.na(rsf_cilla_df$slope),]
Fit a downweighted Poisson regression
m_rsf_cilla <- glm(case_*k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev + slope + var_NDVI,
family = poisson(), data= rsf_cilla_df, weights = w_)
## Warning: glm.fit: fitted rates numerically 0 occurred
Summary of model and confidence intervals for parameter estimates
summary(m_rsf_cilla)
##
## Call:
## glm(formula = case_ * k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev +
## slope + var_NDVI, family = poisson(), data = rsf_cilla_df,
## weights = w_)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.06440 -0.00419 -0.00044 -0.00005 0.59143
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.32153 3.56120 -2.898 0.00375 **
## x_ 0.31218 0.36103 0.865 0.38722
## y_ 0.48704 0.44872 1.085 0.27775
## I(-(x_^2 + y_^2)/2) 1.14409 0.27988 4.088 4.35e-05 ***
## elev -0.01275 0.01742 -0.732 0.46426
## slope -0.03507 0.17592 -0.199 0.84201
## var_NDVI -72.71439 58.86148 -1.235 0.21670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1125.8 on 289061 degrees of freedom
## Residual deviance: 1057.4 on 289055 degrees of freedom
## AIC: 1071.4
##
## Number of Fisher Scoring iterations: 24
confint(m_rsf_cilla)
## 2.5 % 97.5 %
## (Intercept) -16.63420393 -2.92139831
## x_ -0.38125884 1.06078566
## y_ -0.33463784 1.41786282
## I(-(x_^2 + y_^2)/2) 0.67729735 1.77733492
## elev -0.04876778 0.01689897
## slope -0.48964633 0.20232405
## var_NDVI -185.72957712 45.85066051
Map of suitability, including the home ranging behaviour
suitability_glm <- predict_rsf(m_rsf_cilla, cilla_akde, crop(buffalo_env, extent(cilla) * 2))
raster::plot(suitability_glm)
Map of suitability, without the home ranging behaviour
suitability_no_avail_glm <- predict_rsf(m_rsf_cilla, cilla_akde, crop(buffalo_env, extent(cilla) * 2), include_avail = FALSE)
raster::plot(suitability_no_avail_glm)